# LOAD PACKAGES
library("ggplot2")
library("plyr")
library("Rcpp")
library("blockmodels")

########################### ########################### ###########################
# LOAD DATA: European Commission
# Q34a - Internal Contacts 
df_ec <- read.csv(file="./df_ec.csv", header=T, row.names=1, sep=",", na.strings="")

################# BLOCKMODELS EUROPEAN COMMISSION ###########################
# preferred library
library("blockmodels")

df_ec_nw <- as.matrix(df_ec[,c(4:14)])

mymodel <- BM_poisson("LBM", df_ec_nw, explore_min=4, explore_max=10)
mymodel$estimate()
which.max(mymodel$ICL)

mymodel$memberships[[9]]$plot()
mymodel$plot_obs_pred(9)

mymodel$model_parameters[[9]]
mymodel$model_parameters[[9]]$lambda <- mymodel$model_parameters[[9]]$lambda[c(5,1,4,3,2),c(1,3,2,4)]

par(mfrow=c(1,1))
mymodel$plot_parameters(9)

########################### ########################### ###########################

# LOAD DATA: Council Secretariat
# Q33a - Internal Contacts 

df_gsc <- read.csv(file="./df_gsc.csv", header=T, row.names=1, sep=",", na.strings="")

df_gsc_nw <- as.matrix(df_gsc[,-c(1:3)]) #keep only network columns for the blockmodels

################# BLOCKMODELS COUNCIL SECRETARIAT ###########################

mymodel2 <- BM_poisson("LBM", df_gsc_nw, explore_min=4, explore_max=10)
mymodel2$estimate()
which.max(mymodel2$ICL)

mymodel2$memberships[[7]]$plot()
mymodel2$plot_obs_pred(7)

mymodel2$model_parameters[[7]]
mymodel2$model_parameters[[7]]$lambda <- mymodel2$model_parameters[[7]]$lambda[c(2,3,1),c(1,4,3,2)]

par(mfrow=c(1,1))
mymodel2$plot_parameters(7)

